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c~| , We present the first lattice calculation of the B-meson binding energy A and of the 



X 



kinetic energy — Ai/2toq of the heavy-quark inside the pseudoscalar B-meson. This cal- 
culation has required the non-perturbative subtraction of the power divergences present 
in matrix elements of the Lagrangian operator hD\h and of the kinetic energy operator 
hD 2 h. The non-perturbative renormalisation of the relevant operators has been imple- 
mented by imposing suitable renormalisation conditions on quark matrix elements, in the 
Landau gauge. Our numerical results have been obtained from several independent nu- 
merical simulations at j3 = 6.0 and 6.2, and using, for the meson correlators, the results 
obtained by the APE group at the same values of (3. Our best estimate, obtained by com- 
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bining results at different values of /3, is A = 190 MeV. For the MS running mass, 
we obtain m^Wtj,) = 4.17 ± 0.06 GeV, in reasonable agreement with previous determina- 
tions. From a subset of 36 configurations, we were only able to establish a loose upper 
bound on the 6-quark kinetic energy in a i?-meson, Ai = (B\hD 2 h\B) / '(2Mb) < 1 GeV 2 . 
This shows that a much larger statistical sample is needed to determine this important 
parameter. 

* On leave of absence from Dip. di Fisica, Universita degli Studi "La Sapienza", Rome, Italy. 



1 Introduction 



The Heavy Quark Effective Theory (HQET) ||-@ has proven to be an extraordinary tool 
for studying heavy flavour physics. In this approach, physical quantities are expanded 
as series in inverse powers of the heavy quark masses. The spin-flavour symmetries, 
appearing in the infinite mass limit, are then used to relate different hadron masses or weak 
amplitudes which control heavy meson and baryon decays @ . For example, in the infinite 
mass limit, the set of six hadronic form factors, which parameterize the matrix elements 
of the flavour changing vector and axial vector current in B — » D, D* semileptonic decays, 
can be reduced to a single universal one: the so called the Isgur-Wise function Spin- 
flavour symmetries, however, are not sufficient to predict all the properties of the weak 
form factors and of other important quantities such as the meson decay constants and 
the velocity dependence of the Isgur-Wise function. Among the quantities that cannot be 
predicted on the basis of the HQET only, there are several parameters which characterize 
the dynamics of strong interactions, such as the heavy quark binding energy, relevant for 
higher order corrections to the semileptonic form factors, and the heavy quark kinetic 
energy, which enters in the predictions of many inclusive decay rates || . 

The lattice formulation of the HQET offers the possibility of a numerical, non-per- 
turbative determination of these quantities from first principles and without free param- 
eters [J. For example, the most important achievement of lattice simulations of the 
HQET has been the computation of the B-meson decay constant in the static limit, 
fg at . In this work, we present the first lattice calculation of the B-meson binding en- 
ergy, A, and of the kinetic energy of the heavy quark in the B-meson —\\/(2mo), where 
Ai = (B\hD 2 h\B)/(2M B ). 

The parameter A denotes the asymptotic value of the difference between the hadron 
and the heavy quark "pole" mass tuq 

A = lim (M H - m Q ) . (1) 

rriQ^co 

It has been recently shown that the pole mass is ambigous due to the presence of infrared 



renormalon singularities [10, 11 1. At lowest order in 1/mg, the infrared renormalon am- 
biguity appearing in the definition of the pole mass is closely related to the ultra-violet 
renormalon singularity present in the matrix elements of the operator hD^h. This singu- 



larity is due to the linear power divergence of hD^h, induced by its mixing with the lower 
dimensional operator hh. In perturbation theory, using dimensional regularization, the 
power divergence is hidden by the absence of an intrinsic scale in the computation. On 
the lattice, because of the hard cut-off, renormalon poles are absent [O]. In this case, the 
linear divergence manifests itself as a power divergence in the inverse lattice spacing 1/a, 
which appears in the mixing coefficient of the operator hh. In ref. jl^] it was stressed 
that these divergences must be subtracted non-perturbatively since factors such as 

1 / [90(a) dg i 

_ex P(~ / ~ A QCD, (2) 



0(9') 

which do not appear in perturbation theory, give non-vanishing contributions as a — > 



(see also refs. [13, 14]). In this sense, power divergences in theories with a hard cut-off and 
renormalon poles in dimensional regularization are closely related. For a more detailed 
discussion see ref. |tl5| . The intrinsic ambiguity of 0(Aqcd)> present in the renormalisation 
of hD^h, implies an ambiguity in the definition of a finite A and hence of rriQ. 
Falk, Neubert and Luke [16] have proposed a different definition of A, 

x -(0\hTD 4 q\M H } 

(0| hTq \Mh) U 

where h (q) is the effective heavy-quark (light-quark) field, T is a Dirac matrix and Mjj 
a meson annihilated (created) by the operator Jr = hTq (jl = qTh). This definition 
contains the same renormalon ambiguities as that in eq. (||). 

In the lattice HQET, the "binding energy" computed in numerical simulations corre- 
sponds to the definition given in eq. (|j). Consider the two-point function 

C(t) = (0| J r (x, t)4(0, 0) |0> = £ (0| h(x, t)Tq(x, t) q(0, 0)Th(0, 0) |0) (4) 

X X 

For sufficiently large Euclidean time t, 

C{t) -> Z 2 exp(-£t) (5) 

where Z is a constant. The definition (^) implies that A = £: indeed £ can be interpreted 
as the difference Mh — niQ where Mh is the mass of the lightest meson which can be 
created by the operator Jjt. It is clear, however, that £ cannot be a "physical" quantity 
because it diverges linearly as a — > 0. This can be checked in one-loop perturbation theory 
and is a consequence of the mixing of the operator hD^h with hh, as mentioned above. 



It has been argued that it is possible to subtract the divergent term by computing the 
coefficient of hh in perturbation theory |l7[ |. Although it is true that with a hard (i.e. 
dimensional) ultraviolet cut-off, such as the lattice spacing or the Pauli-Villars regulator, 
the matrix elements of the bare operators have no renormalon ambiguities, the subtraction 
of the power divergences using perturbation theory reintroduces renormalons [jl^] . In other 
words, the perturbation series for the power-divergent counterterms contain renormalon 
ambiguities, which, as always, manifest themselves as terms which are exponentially small 
in the coupling constant eq. @. Thus the subtraction of power divergences has to 
be performed non-perturbatively if the resulting matrix elements, such as A, are to be 
unambiguous. 

The matrix elements of the kinetic energy operator, Ai, also contain power divergent 
contributions. In this case, the origin of the divergences is the mixing of hD 2 h with the 
operator hD^h, with a coefficient that diverges linearly, and with the scalar density hh, 



with a quadratically divergent coefficient [12]. Ax determines the 1/mg corrections to the 
heavy quark mass and hence enters many theoretical expressions of weak decay factors. 
As in the case of A, the quadratic and linear divergences of Ai must be subtracted non- 
perturbatively. 

The numerical values of A and Ai, presented in this paper, have been obtained by 



using the non-perturbative method proposed in ref. IS]. In that work, it has been shown 
that a non-perturbative renormalisation prescription, which can be implemented in lattice 
simulations, exists such as to avoid simultaneously both power divergences and renormalon 
ambiguities, in matrix elements and coefficient functions separately. In a theory regulated 
by a dimensionful cut-off, it is consistent not to perform the subtractions of the power 
divergent terms at all, but to work with the bare operators and to compute the coefficient 
functions (which will therefore contain powers of the cut-off) in perturbation theory [O]. 
In this case however, the matrix elements in the effective theory are divergent in the ultra- 
violet cut-off and depend on the regularization. Therefore they cannot be interpreted as 
"physical" quantities, in contrast to the approach that we adopt here. 

The linear divergence in A is eliminated by a suitable redefinition of the operator 
hD^h. This definition corresponds to the same normalization condition that one usually 
imposes in perturbation theory. We require that the matrix element of a combination 



of hD\h and hh, hD^h = hD^h + 5mhh, is zero for given external heavy quark states, 
in the Landau gauge: (h(p4 = Q)\hDlh\h(p4 = 0)) = q Contrary to the perturba- 
tive procedure, which reintroduces the renormalon ambiguities in the matrix elements of 
the subtracted operator, the non-perturbative renormalization condition is unambigous, 
though prescription dependent, and independent of the regularization procedure. It is 
also quite natural in the sense that it allows a "physical" definition of A that is finite and 
independent of the ultraviolet cut-off. This procedure can be extended to the operators 
appearing in higher orders of the 1/rriQ expansion. Moreover it allows the matching of 
the operators of the HQET to those in the full theory (QCD) to be performed, via a 
combination of perturbative and non-perturbative calculations, in such a way that the 
Wilson coefficient functions are free of non-perturbative ambiguities at any given order in 
the 1/rriQ expansion. 

We show below that accurate results are obtained for the binding energy to this order. 
Our best estimate is 

A = 190 +5 ° MeV , (6) 

-30 v ' 

where the error has been obtained by combining the statistical and systematic errors, as 



will be discussed below. Our results show that, as expected jl8|, A is indeed independent 
of the ultra-violet cut-off a -1 , within reasonably small statistical and systematic errors. 

In order to remove the power divergences from the kinetic energy operator, we have 
imposed on the relevant operator a renormalisation condition which corresponds to the 
"physical" requirement {h(p = 0)\hD^h\h(p = 0)) = 0, where hD^h is the subtracted 
kinetic energy operator fTsfl . This renormalisation condition, which will be explained in 
detail in the next section, has been used to extract the values of the mixing coefficients of 
the kinetic energy operator with the lower dimensional ones with a small statistical error. 
Unfortunately, after the subtraction of the power divergences, we were unable to obtain 
a precise value for Ai, because of the large cancellations between the operator matrix 
element and its counterterm. We can only put a loose upper bound of 1 GeV 2 on Ai. 
Nevertheless, the results of this study are so encouraging that we are implementing this 
procedure on the APE100 computer to perform a high statistics lattice calculation of both 



1 This requires certain assumptions on the infrared behaviour of the heavy quark propagator that will be 
discussed below, see also fllal. 



A and Ai, whose results will be published elsewhere. We believe that the present results 
demonstrate the feasibility of the method proposed in ref. jl8| to compute quantities 
relevant in heavy flavour phenomenology. In this way, it will still be possible to use the 
HQET and the notion of a "pole" mass, now defined non-perturbatively, which seemed 
to be ruined by the presence of the renormalons. Preliminary results of the present study 
can be found in ref. [pUf! . 

The plan of the paper is the following. In sec. [2| we introduce the relevant formulae 
which define the non-perturbative procedure for renormalising the operators hD^h and 



hD 2 h [18]; in sec. |we describe the numerical calculation of A and Ai and discuss the main 
results of this study; in the conclusion we present the outlook for future developments 
and applications of the method discussed in this paper. 



2 Non-perturbative definition of A and A 



1 



In this section we define the renormalisation prescription which we will use to calculate 
"physical" values of A and Ai. The prescription involves imposing appropriate renormal- 
isation conditions on the quark matrix elements of the operators hD^h and hD 2 h, such 
that all their matrix elements are free of power divergences flq]. Similar methods have 



been used for light quark operators in refs. 2C| 



In numerical simulations, quark and gluon propagators can be computed non-pertur- 



batively by working in a fixed gauge, typically the Landau gauge |20|-|23|]. The heavy 
quark propagator, at lowest order in 1/toq, has the form 

S(x,t) = {S(x,t\6,0)) = 5(x)9(t)5 ij A(t) exp(-At) (7) 

where i,j are colour indices; 

S(x,t\y,w) = 5(x -y)9{t - w) expfi / A (t')dt') (8) 

is the non-translationally invariant propagator for a given gauge field configuration, com- 
puted in a given smooth gauge, typically the Landau gauge, and (. . .} represents the 
average over the gauge field configurations. A(t) is an unknown function of t, and we 
assume that it decreases more slowly than an exponential at large times, specifically we 



require that 



lim 



a V 



A(t + a)\ d , . 

^ M ~ lim — ln^(i) = 0. 
A{t) ' t->oo dt y ' 



0) 



Below we will show that, within the precision of our simulations, our results for the heavy 
quark propagator are consistent with the condition in eq. @. As will also be explained 
below, the condition @ is not strictly required for the definition and determination of A, 
since this can be done using the values of the propagator at small times t. However in 
that case it no longer has a direct interpretation as a binding energy, and our preferred 
definition of A does use the behaviour of the propagator at large t. 

The constant A in eq. (0) is linearly divergent in 1/a and would correspond, in dimen- 
sional regularisation, to an ultraviolet renormalon in the effective heavy-quark propagator. 
Since the linear divergence in A is due to the mixing of the operator hD^h with the con- 
served scalar density operator hh, we can remove it by adding to the Lagrangian of the 
lattice HQET 

jr cS = h( x )D A h{x) (10) 
a counter-term of the form 5mh(x) h(x). The HQET Lagrangian then becomes 

4 ff = —=—(h(x)D 4 h(x) + 5mh(x)h(x)), (11) 

1 + a dm v / 

where the factor 1/(1 + a 5m) has been introduced to ensure the correct normalization of 
the heavy quark field h. With the action C' eS , the heavy-quark propagator is given by: 

ln(l + a 5m)' 



S'(x,t) = 5{x)6{t)5 l] A{t) exp 



A + 



(12) 



The mass counter-term is defined by the behaviour of S(x, t) at large values of the time 



ln(l + a 5m) 1 

5m = = lim — in 

a t^co a 



Tr (S(x, t + a 
Tr(s(a?,f)) 



lim 

t— >oo 



a V A(t) 



A 



-A, 



(13) 



where the traces are over the colour quantum numbers, and we have assumed the validity 
of the condition in eq. (^). Our numerical results, support the validity of this condition 
and the use of eq. (Il3|) is our preferred determination of 5m. 



We now define the renormalised binding energy by 

A = S - 5m, (14) 

which corresponds to the following relation between the meson and the heavy quark mass 

M H = niQ + £ - 5m (15) 

rriQ can be interpreted as a subtracted pole mass, and contains no renormalon effects. A 
similar relation can be found in the case of a heavy baryon. 

The definition of 5m given in eq. (|i~3| ) is not unique. A possible alternative definition 
would be, for example, 



5m(t*) = — In 
a 



Tr 



1 , [A(t* + a)\ 



Tr(s(x, t*) 

where t* is a given time at which we perform the subtraction. The corresponding defi 



nition of A, see eq. (14) above, will clearly depend on the choice of t*: t* parametrizes 
the renormalisation prescription dependence and can be considered as the renormalisa- 
tion point in coordinate space. For physical matrix elements, the residual mass appears 
only through the combination ttjq — 5m, in such a way that different choices of 5m are 
compensated by different values of mg The use of the propagator at small times, 

t*A.QCD *C 1, to define 5m(t*), and hence A(t*) does not require any assumption about 



the behaviour of A(t) at large times, and in section 3.2 we present the results for A(t*) 
obtained in this way. 

In addition to the non-perturbative contribution of O(Aqqd) to 5m(t*), there is a 



perturbative one proportional to 1/t* [18], 

-^ pcrt (O = -^^ + 0(^) (17) 

where 7^ is the one-loop contribution to the anomalous dimension of the heavy quark 
field (70 = —6 in the Landau gauge) and Cf is the quadratic Casimir operator in the 
fundamental representation (Cf = 4/3). Thus the definition of A(t*) defined, at small 
times cannot readily be identified as a physical binding energy. Nevertheless, computed 
values of A(t*) can be used to determine standard short-distance heavy quark masses (such 



as the MS one) using perturbation theory (as will be explained in section 3.1 below). This, 
together with the fact that no assumption about the infra-red behaviour of the heavy quark 
propagator is necessary, is of fundamental importance. 



2.1 Non-perturbative subtractions for Aj 

The renormalised kinetic operator hDgh, free of power divergences, has the form 

hix) 31 hix) = hix) 3 2 hix) - — - — -— r- f hix) D 4 hix) + 5m hix) hix) ) 

a ( 1 + a dm ) V / 



^h(x)h(x), (18) 



a (1 + a 5m) 

h It I h ( nr 1 

I 2 

where the constants ci and C2 are functions of the bare lattice coupling constant go{a). 
They have been computed in one loop perturbation theory in ref. [^]. Notice that we 
have preferred to express hD$h in terms of the subtracted operator which explicitly 
contains the residual mass 5m. In this way we can use the equations of motion of the 
Lagrangian C! given in eq. (]TT|). This will prove useful below. The constant C2 enters in 
the renormalisation of the heavy quark mass. Therefore, it will contribute to the relation 
between Mh, mQ and A at order 1/mg (see below). On the other hand, the constant 
ci contributes to the renormalisation of the heavy quark wave-function and hence to the 
renormalisation of all the operators containing a heavy quark field, but not to the relation 
for the quark mass. 

In order to eliminate the quadratic and linear power divergences, a possible non- 
perturbative renormalisation condition for hDgh is that its subtracted matrix element, 
computed for a quark at rest in the Landau gauge, vanishes 

(h(p = 0)\hD 2 s h\h(p = 0)) = 0. (19) 

This is equivalent to defining the subtraction constants through the relation (in the fol- 
lowing we will work in lattice units, setting a = 1) 

Pniit) = c i + c 2t, (20) 

where 

_ EAS a '(x,t\o,Q)) = Eko Ess j s '& % g) & y (t' ) SU f\6, o) ) 
P& Ex(s>(x,t\d,o)) EAS'(x,t\o,o)} 

By fitting the time dependence of pg 2 (t) to eq. (|20|), one obtains c\ and c^. 

The heavy-quark propagator that enters in eq. ( |2l| ) is the subtracted one, i.e. it is 
calculated with the action (11) instead of (0). We now demonstrate that p^it) can be 
expressed in terms of unsubtracted propagators only: 



Pd^) 



Eko Ex,y ((S'(S, t\y, t>) D 2 y (t>) S'(y, t% 0))) 
Z x {S>(x,t\6,0)) 



Eko E^((5(x,i|y,O^ ( *- t ' ) ^(<0^(y ) i / |0,0)e^ t ')) _ 

(22) 



t|ff,0) e Smt ) 
n=o Eg,y ((S(x, t\y, t>) Dffl) S(y, t'\0, 0))) 



Eg(S(x,t\o,o)) 

Notice that this argument holds for any operator which does not contain a time derivative. 

For some important applications it is only the constant C2 which is required. C2 can 
also be determined directly by eliminating the sum over t' in eq. (f2~l|): 



M H = m Q + S-Sm-(l-^X m )[^- -* +0(- r ), (24) 



, +l T,x,y ((S'(S, t\y, t') D 2 y (t>) S'(y, t'\0, 0))) 

C2 = p »™ = EAS^Mo)) (23) 

for t' + 0, t. 

The relation between the mass of the meson and the mass of the quark to order 1/mg 
is then given by 

4^1 J ^ 2m Q J + X 
where A^ 6 = {B\hD 2 h\B)/(2M B ). \\ am can be determined from a computation of two- 
and three-point correlation functions in the standard way. Consider the meson three-point 
correlation function (the extension of this discussion to baryons is entirely straightforward) 

C B2 (t',t) = J^(0\J T (x,t) h(y,t')D 2 y h(y,t')4(0,0)\0) (25) 

For sufficiently large values of t' and t — t' 

Cg 2 (t', t) -» Z 2 \\ &ve exp (-(£ - 5m)i) . (26) 

A convenient way to extract A^ are is to consider the ratio 

^',*) = %!^-A^ (27) 

As usual Xi are must be evaluated in an interval in which R(t', t) is independent of the 
times t' and t, so that the contribution from excited states and contact terms can be 
neglected. 

The term proportional to Xg 2 in eq. (24) is absent in continuum formulations of the 



HQET, and is a manifestation of the lack of reparametrisation invariance in the lattice 
version. It has been calculated in ref. | |12| |. Notice that only the constant C2 enters the 
eq. ( p4| ) because c\ is eliminated by using the equations of motions. c\ only modifies 
the wave function renormalisation of the heavy quark, thus contributing to the 0(l/mg) 
corrections of the hadronic matrix elements. 



simulation 


volume 





Number of configurations 


set A 


16 3 x 32 


6.0 


36 


set B 


16 3 x 32 


6.0 


300 


set C 


20 3 x 32 


6.2 


50 


set D 


18 3 x 64 


6.0 


210 


set E 


18 3 x 64 


6.2 


420 



Table 1: Parameters of the numerical simulations, the results of which have been used for 
present study. 

3 Numerical implementation of the renormalisa- 
tion procedure 

As explained in the previous section, the determination of A and Ai requires the compu- 
tation of the quark propagator and matrix elements between quark states in a fixed gauge 
(in order to obtain the subtracted operators), as well as the evaluation of matrix elements 
between hadronic states. We have obtained our results using five independent numerical 
simulations, whose main parameters are given in table [[]. 

Our best value of the subtracted binding energy, A = £ — 5m, has been determined 
by combining the values of 5m obtained using set B and set C with the calculation of £ 
performed by the APE collaboration at j3 = 6.0, set D , and 6.2, set E [2(|. £ had been 
determined using the SW-Clover fermion action for the light quarks. The calculations were 
performed at several masses of the light quark, so that extrapolations to the chiral limit 
are possible. We also present the results for the subtraction constants c\ and C2 obtained 
with set B and set C . 

So far we have only computed A^ 3 * 6 using set A . Again, for the light quarks the im- 



proved SW-Clover action [27] was used in the quenched approximation. These exploratory 
calculations were performed at one value of the mass of the light quark, k = 0.1425, for 
which the mass of the corresponding "pion" is about 900 MeV. The details of the simu- 
lation can be found in refs. Preliminary results for both A and Ai for mesons, 
evaluated using this dataset have been presented in ref. [[H]] . 

All the errors have been computed with the jacknife method by decimating one con- 



figuration at a time ( set A and set C ) or five configurations at a time ( set B ). The error 
on £ was computed with the jacknife method also and we refer the reader to refs. (25|, E(| 
for details. 



3.1 Determination of the residual mass 5m 

A possible lattice expression for the forward heavy-quark propagator, to leading order in 
the heavy quark mass, is given by 

S(x\0) = <5(f)#(x 4 )7^(x 4 |0) (28) 

where Vg (x 4 | y 4 ) is the lattice path ordered exponential from (x, y 4 ) to (x, x 4 ), cf. eq. (||), 
usually called "P-line" , 



VA^\y l ) = n u ^ g - x " - n °). * A> y' 

n=l 

V x {x A \y A ) = 1 x A = y A (29) 



This propagator corresponds to the following choice for the covariant time derivative, 
D 4 f(t) = l/a (f(t) - U\(t - a)f(t - a) ). 

In order to reduce the statistical noise, we have computed, in the lattice Landau gauge, 
the quantity 

S H(t) = ^E(^[^ 4 = t l°)])' ( 3 °) 

X 

where the trace is over the colour indices and V denotes the spatial volume of our lattice. 
It is this averaged propagator Suit), which has been used in the computations below. 

There is a subtle point that we would like to discuss briefly. It can be demonstrated 
that 0(a) effects in heavy-light operator matrix elements between physical states are 
cancelled by improving the light quark propagators only |28|]. On the other hand, in order 
to improve off-shell matrix elements, which is the case when renormalising the operators 
between quark states, it is necessary to use an improved version of the heavy quark 
propagator in the effective theory, for example 

•IN y 4 -IE 4 + l' 



1 



7>,(x 4 |y 4 ). (31) 



Notice that the improved P-line tends very rapidly to the unimproved one as x 4 — y 4 
increases. The propagator in eq.(|3l|) corresponds to the following time derivative Z?4/(i) = 
1/a (3/2/(t) - 2£/t(i - a)f(t - a) + l/2U*(t - a)U\t - 2a)/ (t - 2a) ). It is also possible 
to add a residual mass term to the heavy quark action in such a way that it modifies the 
propagator ( |3~l| ) by an exponetial in time (up to an overall normalisation factor). Such as 
mass term takes the form 3/2 (1/ X-l)h(t)h(t) + l/2 (A- l)h(t)h(t-2a). In the following, 
when discussing the improved heavy quark propagator, we will implicitly assume that the 
mass term is of this form. 

To determine the residual mass, we have to compute the effective mass of the propa- 
gator Sn(t), defined by 

„ M() ._ ln (M^)) (32) 

In figs. [I] (from set B ) and |2] (from set C ), we present the values of 5m(t) for the improved 
and unimproved propagators as a function of t/a. The effective mass is indistinguishable 
in the two (improved and unimproved) cases, for t/a > 4-5. Thus, in order to minimize 
lattice artefacts, we have only used the results obtained for t/a > 5. Inspired by the 
results of one-loop perturbation theory we made a fit to 8m{i) using the expression 

a 6m(t) = a 8m + 7— (33) 

In order to mimic higher order effects, we have also used different expressions to fit 5m(t), 
e.g. 

a 5m(t) = a 5m + 7' ln^— - — ^ (34) 



or 



ff1 ( a s [1C/{t + a)} ^ „ flu[{t + a)]+C^ 

a5m{t) = a5 m - 1 ln( ^ )^a5m + 1 ln( ^ + c ), (35) 

and changed the interval of the fits in order to check the stability of the determination 
of 8m. In eqs. (p3|)-(35), 8m, 7, ... ,7" and C are free parameters of the fit. The curves 



shown in fig. || and ^ correspond to fits of the improved heavy quark propagator to eq. 
(H), in the interval 5 < t/a < 12. 

From the different results obtained by varying the fitting functions and the time in- 
tervals, see tables |2| and ||, we quote 

a8m = 0.521 ± 0.006 ± 0.010 at /3 = 6.0 (36) 



Figure 1: Effective mass of the heavy-quark propagator Suit), at {3 = 6.0, as a function of the 
time. The curve represents a fit of the numerical results (in the improved case) to the expression 
given in eq. (34). 



Figure 2: Effective mass of the heavy-quark propagator Sjj{t), at (3 = 6.2, as a function of the 
time. The curve represents a fit of the numerical results (in the improved case) to the expression 
given in eq. $Jd]). 



Subtraction constant aSm at (3 = 6.0 


Fit 


t = 4 - 12 


t = 5 - 12 


t = 5 - 14 


t = 6 - 14 


t = 7- 14 


t = 8 - 14 


eq. (33) 
X 2 /dof 


0.507(4) 
1.50 


0.515(6) 
0.56 


0.515(6) 
0.86 


0.52(1) 
0.84 


0.52(2) 
0.96 


0.51(3) 
1.05 


eq. (34) 
X 2 /dof 


0.515(4) 
0.97 


0.521(6) 
0.46 


0.521(6) 
0.81 


0.53(1) 
0.85 


0.53(2) 
0.98 


0.51(3) 
1.05 


eq. (35) 
x7dof 


0.513(4) 
1.20 


0.521(6) 
0.56 


0.520(6) 
0.93 


0.51(3) 
0.97 


0.51(1) 
1.10 


0.50(2) 
1.28 



Table 2: Numerical values of the constant a5m found by using the results of set B , at (3 = 6.0. 
The results are from several fits in different time intervals. We also give the uncorrelated 
X 2 /dof in the different cases. The numbers given in this table refer to the improved heavy quark 
propagator only. 



Subtraction constant a5m at (3 = 6.2 


Fit 


t = 4 - 12 


t = 5 - 12 


t = 5 - 14 


t = 6 - 14 


t = 7- 14 


t = 8 - 14 


eq. (33) 
X 2 /dof 


0.437(5) 
0.50 


0.441(7) 
0.93 


0.440(8) 
0.80 


0.45(1) 
0.75 


0.46(2) 
0.70 


0.45(2) 
0.83 


eq. (34) 
x7dof 


0.442(5) 
0.40 


0.445(8) 
0.85 


0.445(8) 
0.74 


0.45(1) 
0.73 


0.46(2) 
0.70 


0.46(3) 
0.83 


eq. (35) 
x7dof 


0.443(6) 
0.50 


0.445(8) 
1.02 


0.445(8) 
0.85 


0.45(1) 
0.85 


0.45(5) 
0.84 


0.46(3) 
1.05 



Table 3: Numerical values of the constant aSm found by using the results of set C , at (3 = 6.2. 
The results are from several fits in different time intervals. We also give the uncorrelated 
X 2 /dof in the different cases. The numbers given in this table refer to the improved heavy quark 
propagator only. 



a5m = 0.445 ± 0.008 ± 0.010 at = 6.2 (37) 

where in both cases the first error is statistical, and the second is an estimate of the sys- 
tematic uncertainty, based on the spread of results obtained using different time intervals 
and fitting functions. 

The determination of the mass counter-term at fixed t = t*, requires no fitting, and 
the results obtained using set B and set C are presented in table ||] below. 



3.2 Determination of A 

We are now ready to present our prediction for A. In order to evaluate the subtracted A, 



we have used the results of the high statistics calculations of £ given in refs. [25, |26| ( set D 
and set E ). We will also make use of the results obtained by using the standard Wilson 
action, on a 18 3 x 64 lattice, at = 6.0, with a statistical sample of 200 configurations p5|]. 
In order to obtain A we have used: 



• 5m from eqs. (|36|) and (|37|); 

• the SW-Clover determination of £ of the APE collaboration, a£ = 0.61 ± 0.01 at 



= 6.0 and a£ = 0.52 ± 0.01 at = 6.2 |25], |§; 

• a~ x {0 = 6.0) = 2.0 ± 0.2 GeV and a~ 1 (0 = 6.2) = 2.9 ± 0.3 GeV. The calibration of 
the lattice spacing in quenched simulations typically has an uncertainty of O(10%), 
depending on the physical quantity which is used to set the scale. We take these 
results as a fair representation of the spread of possible values. 

We then find 

A = £ - 5m = 180 ± 35 MeV at = 6.0 (38) 
A = £ - 5m = 220 ± 55 MeV at = 6.2 (39) 

where the statistical errors have been combined in quadrature with those due to the 
uncertainty in the lattice spacing. 

Within the uncertainties, the results in eqs. (|38|) and p9| ) are compatible with the 
expected independence of A of the lattice spacing. Given the intrinsic uncertainty in the 
value of the lattice spacing in quenched simulations it is difficult however, to check this 
more precisely. Indeed, using different physical quantities to set the scale can increase 



or decrease the difference in the central values of A at (3 = 6.0 and 6.2. For example, 
using the string tension to set the scale one finds a~ l ((3 = 6.0) = 1.88 GeV and a~ l ((3 = 
6.2) = 2.55 GeV, giving A = 170 ± 30 MeV at (3 = 6.0 and A = 190 ± 40 MeV at 
,9=6.2, whereas using the mass of the p-meson to set the scale the APE collaboration 
finds a" 1 (/3 = 6.0) = 1.95 ± 0.07 GeV and a - l ((3 = 6.2) = 3.05 ± 0.20 GeV [§5], ||], 
which corresponds to A = 176 ± 30 MeV at (3 = 6.0 and A = 228 ± 50 MeV at /3 = 6.2 |. 
Nevertheless in both cases the results are compatible at the two values of (3. 



Assuming that A is indeed constant in a, we combine the results in eqs. (^) and (3£) 
to obtain 

A = (190 ± 30) MeV . (40) 

Before quoting our final result, we need to estimate the discretisation error. 

From a comparison of the values of 5m obtained with the improved and unimproved 
heavy quark propagators, we believe that discretisation effects are negligible for this quan- 
tity. Indeed discretisation errors in quantities which only depend on the gauge fields are 
of 0(a 2 Aq CD ), when evaluated using the Wilson gauge action. However, in the compu- 
tation of A (and Ai), correlation functions which contain the light quark propagator are 
evaluated, and with the SW-Clover and Wilson fermion actions this introduces errors of 
0(a s a Aqcd) and 0(a Aqcd) respectively. Notice that these effects are formally larger 
than the higher order corrections to A (we work in the approximation a -1 <C itiq). 

To obtain an estimate of the size of the discretisation errors, we compare £ obtained with 
the standard Wilson action and the SW-Clover action at the same value of (3, (3 = 6.0. 
In the Wilson case, by working at four different values of the light quark mass, the bare 
binding energy, extrapolated to the chiral limit in the light quark mass, was found to be 
a£\y = 0.608(8). In the SW-Clover case, by working at three different values of the light 
quark mass, the result for the bare binding energy, extrapolated to the chiral limit in the 
light quark mass, was found to be a£$w = 0.616(4) The difference between the central 
values obtained with the two actions £sw — &w = (0.616 — 0.608) a -1 ~ 16 MeV. We 



2 Notice that, using the mass of the p-meson, the UKQCD collaboration found a 1 (/3 = 6.2) ~ 2.7(1) GeV pQ |, 

corresponding to A = 203 ± 45 MeV. 

3 As a check of our calculations with set A . we have verified that on these configurations, at the value of the 

mass where we have computed the light quark propagator (K = 0.1425), our results for £ agree with those of 

refs. 



deduce that +20 MeV is a reasonable estimate of the discretisation error in the determi- 
nation of A. We therefore quote as our final result for A 

A = 190 +5 ° MeV . (41) 

-30 v ' 

The prediction given in eq. (|4ll ) can be compared with other results that have been 
presented in the literature. In perturbation theory one finds 

ot f d^o / 1 \ 

a6m ^ = f J w (gu^w = 2-12 x as (42) 

By using values for a s which are commonly proposed in the literature for the "boosted" 
coupling ||, 0, a s = 0.13-0.18 at p = 6.0, eq. @ would give a5m pcrt = 0.28 - 0.38. 
Thus, even in the most favourable case, Jm pert is about 280 MeV smaller (i.e. A per t is 
about 280 MeV larger) than our non-perturbative determination. 

In ref. |3l[], the bare binding energy £ has been determined, using the Wilson action 
for light quarks, on a variety of lattice volumes and at several values of /3, f3 = 5.7, 5.9, 6.1 
and 6.3. The results are consistent with a linear dependence 

a£{a) = £ + a A FNA l (43) 

where £q and Afnal are parameters of the fit, £q = 0.351(14) and Afnal = 0.481(25) GeV. 
The value of £q is consistent with a5m pert computed using an "effective" a s = 0.166 (this 
value may be considered as an average of the values of the strong coupling constant on 
the points in (3 where £ has been computed). On the other hand, the value of the "fi- 
nite" binding energy Afnal is about 300 MeV larger than ours 0. Our interpretation is 
that, up to possible 0(a) effects, the two determinations differ because of the finite non- 
perturbative contribution of O(Aqcd) that has been subtracted only in our case. Using 
the definition of ref. [31| however, it is not clear how to match the full and the effective 



theories, since their definition includes non-perturbative, uncalculable effects. 

A further demonstration of the existence of the non-perturbative effects is provided by 
the comparison of a5m and a -1 at /3 = 6.0 and 6.2. In the absence of non-perturbative 
terms of 0(Aqcd) 5 i-e. if 5m is given only by the linearly divergent contribution, we should 
find R m = a5m(p = 6.Q) /aSm(fi = 6.2) ~ R as = a s {(5 = 6.0) /a a (J3 = 6.2). Numerically 



4 Given the presence of terms of 0(a) the stability of the results with respect to a quadratic fit of the form 
a£(a) = £ o + Afnal^ + ^2a 2 , where £% is a constant, remains to be checked. 



R 
P 


i /a 


CLomyz ) 


7T/ I [KjG V J 


A{t ) ^iviev j 


6.0 


3 


0.3670(6) 


2.1 ± 0.2 


490 ± 20 ± 50 




4 


0.3980(8) 


1.6 ± U.2 


/ion i on i /in 

420 ± 20 ± 40 




5 


0.4177(9) 


1.3 ±0.1 


390 ± 20 ± 40 




6 


0.4328(13) 


1.0 ±0.1 


350 ± 20 ± 40 


6.2 


4 


0.3484(13) 


2.3 ±0.2 


500 ± 30 ± 50 




5 


0.3663(16) 


1.8 ±0.2 


450 ± 30 ± 50 




6 


0.3773(20) 


1.5 ±0.2 


410 ± 30 ± 40 




7 


0.3842(27) 


1.3 ±0.1 


390 ± 30 ± 40 



Table 4: Results for A(t*) = £ — 5m(t*) for different normalization times t* , using the results 
from set B-set E . The first error on A(t*) is obtained by combining the errors on 8 and 5m(t*) 
in quadrature; the second error (and the error on ir/t*) comes from the calibration of the lattice 
spacing. 



we find R m = 1.16(4) to be compared with R as = 1.03 — 1.06: 1.03 is simply 6.2/6.0; 1.06 
has been estimated from R a es, where a| ff = a^ att / (□) with a^ att = (6//3)/(47r) and (□) is 
the expectation value of the plaquette. 

We now present the results for A defined at a fixed value of t* (A(i*)). In order to 
be able to use perturbation theory to determine values corresponding to standard short 
distance definitions of the heavy quark mass, t* must be chosen to be sufficiently small. 
In table Q we present the results for the mass counterterm 8m(t*) at both (3 = 6.0 and 
6.2, obtained using the configurations of set B and set C respectively, for small values 
of t* . We then combine the results for 5m{t*) with those for S obtained by the APE 
collaboration ( set D and set E ) to obtain A(£*). These results for A(t*) will be used in 



section [3^ to determine the MS mass. 

We end this subsection with an obvious but important remark. A can be defined in 
many different ways, which correspond to different renormalization prescriptions for the 
renormalized hD^h operators. The presentation of results or bounds for A must therefore 
be accompanied by the definition of the prescription to which they correspond. 



3.3 The MS mass of the 6-quark 

We now give the relevant formulae necessary to match the subtracted mass of the quark 
Wg to the running mass niQ, computed in the MS scheme at the scale \x = mg. We 
introduce the following quantities 

mg(t*) = M H - A(t*) = M H - £ + Sm(t*) , (44) 

C m (f) = 1 - _ _L Urnit*) _ 0,(0)-) . (45) 

37r mQ V a J 

It is straightforward to show, that, at order a s , the relation between Tuq and m,g(t*) is 
given by 

raQ = ra^(f)xC m (f). (46) 
Equation (E^) holds also for t* — > oo, provided at the same time <5m(i*) — » 5m. 



At this order in a s we can write 

»«g = (^-£ + a.(a)f)(l-^9)) (47) 
which is equivalent to the procedure where the linearly divergent term is subtracted in 



perturbation theory. Now both factors in eq. (47) contain renormalon singularities. We 



can also rewrite eq. (47) in the form corresponding to the procedure where we use the 
unsubtracted, linearly divergent "pole" mass Mjj — £, 

m Q = (M H -e)(l- + a a {a)JL) , (48) 

^ ' ^ 6tt arriQ' 

where the divergent dependence on a in the pole mass is compensated by that in the 
coefficient function. In this case no renormalon singularities arise in higher orders, but 
the unsubtracted pole mass and the coefficient function both contain power divergences. 
As required, the relation (|48| ) is independent of the subtraction constant 8m(t*). 

In the numerical evaluation of mg from eqs. ( |4"6| ) and ([47]), we used Mb = 5.278; 



£ from the APE results, see subsection 3^; a5m(t*) from eqs. ([?(]) and ( p7| ) and table 
a s (a) was taken in the range 0.13 and 0.18. To obtain a distribution of values, we 
varied £, Aqcd and aSm(t*) according to a gaussian distribution; a -1 was varied with 
flat distribution within its error, while a s (a) was written in terms of the leading quenched 
expression of the running coupling constant, evaluated at the scale ir/a, with Aqcd dis- 
tributed according to a flat distribution of width a and such that a s = 0.13 for Aqcd — cr 



Figure 3: The ditribution of values of the MS b-quark mass mi,, for t* — > oo ; at (3 = 6.0 and 
(3 = 6.2. Similar distributions are obtained for m b using t* = 3-7. 

and a s = 0.18 for Aqcd + f- The resulting distribution is a pseudo-gaussian, as can been 
seen from fig. || where two histograms of values of m&, corresponding at /3 = 6.0 and 6.2, 
are shown. From the width of the distribution we estimate the average value and error 
on mj. Using eq. fl4q) , we obtain at (3 = 6.0, 

m b = 4.18 ± 0.07 GeV for i* -> oo (49) 

and 

m b = 4.21 ± 0.07 GeV for t* = 3 - 6 . (50) 
The corresponding numbers at f3 = 6.2 are 

m b = 4.11 ± 0.09 GeV for t* -> oo (51) 

and 

m b = 4.13 ±0.09 GeV for t* = 4- 7 . (52) 



Using eq. ([17]), we obtain instead 

m b = 4.22 ± 0.07 GeV at p = 6.0 , 

m b = 4.15 ±0.08 GeV at /? = 6.2 (53) 

The difference between the results of eqs. (^9|)-(|52]) and ( |53"| ) can be interpreted as due 
to higher order corrections in a s . By combining the above results together we estimate 

mb = 4.17 ± 0.05 ± 0.03 GeV (54) 

where the second error is the systematic error coming from the different methods used to 
extract rn b at this order in a s . 



3.4 Determination of c\ and C2 

We have computed the ratio pg 2 (t), defined in eq. (^), using unsubtracted heavy-quark 
propagators, as explained in section 2.2. In order to do this calculation, we need the 
expression of the heavy-quark propagator with the insertion of hD 2 h 

S a (x\y)= J2 S(x,\w) Dl(w 4 )S(w\y) , (55) 



where the lattice heavy quark propagator S(x\w) has been defined in eq. (p8[), and in 
the improved case, we have used the definition of the P-line given in eq. (|3l|). For the 
discretised version of D 2 we have taken 

( 3 

...( 



Dl 



4 E (U k a(3 (x)6 xx + ak + U k J(x - ak)6 xx _ ah - 26 a ,8 x , x ) , (56) 
a k=i 



In fig. we plot pg 2 (t), as defined in eq. (21), as a function of the time t, at /3 = 6.0 
from set B . In the same figure, we also give the result of a linear fit of pgiif) to eq. (2C) 
in the interval 6 < t/a < 12. Similar results were obtained at /3 = 6.2 using the data of 
set C . The dependence of pg 2 (t) on t is in remarkable agreement with the predicted linear 
behaviour Q In order to monitor the stability of the results, we have fitted pg 2 (t) using 
different time intervals and in table || we show our results for c\ and C2, in the improved 



5 As expected, we found that the results for C2 in the improved and unimproved case are completely 
patible. The latter are not reported here. 






time interval 


Cl 


C2 


X /dof 


6.0 


4-12 


0.06(2) 


-0.759(6) 


1.24 


6.0 


5-12 


0.01(5) 


-0.748(10) 


0.98 


6.0 


5-14 


0.01(5) 


-0.748(10) 


1.09 


6.0 


6-14 


-0.16(13) 


-0.724(22) 


0.80 


6.2 


4-12 


0.10(4) 


-0.698(9) 


1.01 


6.2 


5-12 


0.16(9) 


-0.708(18) 


0.98 


6.2 


5-14 


0.17(9) 


-0.710(17) 


0.95 


6.2 


6-14 


0.37(18) 


-0.739(29) 


0.62 



Table 5: Results for the improved renormalisation constants of the operator hD 2 h obtained by 
a linear fit to pg 2 (t) . The time interval of the fit is also given. 



Figure 4: The ratio pg2{t) (improved case) as a function of the time, at (3 = 6.0 from set B . A 
linear fit of pg 2 {t) to the expression in eq. $2(\) in the interval 6 < t/a < 12, is also given. 



Figure 5: The ratio pg 2 (t',t) defined in eq. (p?5|j, in the improved case, at (3 = 6.0 from set B , 
as a function of t' , at several values oft, t = 6-10. 



case, at (3 = 6.0 and 6.2. At = 6.0, we observe a shift of the value of C2 towards 
smaller values as we increase the minimum t-distance (t m i n = 4, 5, 6) at which the fit is 
performed. Since at /3 = 6.2, we find the opposite behaviour, i.e. the value of C2 is shifted 
towards larger values as t m i n is increased, we believe that the shift is a statistical effect 
rather than a systematic one. From table [5], we also observe that it is very difficult to 
determine the value of c±, which, for the improved propagator, seems to be small, with a 
large relative error, and is very unstable with respect to a change of the fitting interval. 
We expect that this instability, which is correlated to the shift of the value of C2 with 
tmim will be reduced with more accurate data for pg 2 (t)- Notice that c\, unlike £ and C2 
which are long-distance quantities, depends on the lattice regularization, i.e. is different 
for the unimproved or the improved heavy quark propagator. 

We also present the results for C2 = p^ 2 {t f ,t), at j3 = 6.0, obtained by using eq. (|23|) . 
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3-5 
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-0.693(18) 


0.40 


6.2 


10 


3-7 


-0.723(29) 


0.22 



Table 6: Results for the renormalisation constant computed from a weighted average 
Pfrit', t) in t' , at fixed t. 

In fig. ||, we show p^ 2 (t',t), as a function of < t' < t, at several fixed values of t, 
t = 6-10. Up to contact terms, we expect pg 2 (t',t) to be a constant in t', at fixed t, 
and also to be independent of t. If the contact terms were entirely due to the mixing 
of the kinetic energy operator with the inverse propagator, eq. (|l~8|), we should find two 
spikes, at t' = and t' = t, and a constant value of pg 2 (t',i) for t' ^ 0,t. The presence 
of operators of higher dimension, due to discretisation errors, introduces terms which 
behaves as derivatives of 6- functions (in time), giving rise to the bell-shape behaviour of 
pg 2 (t',t) shown in fig. [|. Thus in order to obtain C2, we have to look for a plateau in 
the central region in t' , at large values of t. From the figure, we see that it is possible to 
recognize a plateau in t' for t = 8-10. At values of t smaller than t = 8, the contact terms 
are visible at all values of t'; at values of t larger than t = 10 the statistical error become 
quite large. There is a slight shift towards larger values of C2 as t is increased. As discussed 
above, since the effect is opposite at 3 = 6.2, we do not believe that this is a systematic 
effect. In table |6| we present the values of C2, computed from a weighted average of values 
of pg 2 (£'> t) for different t' , at fixed t. The average has been performed only in the central 
region, where there appears to be a plateau. For the sake of comparison, we present the 
results for several values of t, including small ones. 



From the results given in tables [| and ||, and taking into account the previous discus- 
sion, we believe that the best estimate of C2 is obtained from pg 2 (t',t), with t = 8 and 
f = 3-5 

c 2 = -0.73 ±0.01 ±0.02 at [3 = 6.0, (57) 

c 2 = -0.68 ±0.01 ±0.02 at [3 = 6.2, (58) 

where the second error comes from the variation of the values of C2 with t. These results 
can be compared with perturbation theory, which gives C2 = —5.19 x a s ~ — (0.67- 
0.93) for a s = —.13-0.18. For c\ using the improved propagator, such a comparison is 
impossible, due to the relatively large uncertainties in the non-perturbative determination. 



3.5 Determination of the kinetic energy Ai 

The results for Ai have been obtained with limited statistics, using the data of set A . As 



explained in subsection 2A, the value of Xi are can be obtained from R(t',t) as defined in 
eq. (p7|). In principle, we should evaluate R(t',t) using the subtracted propagators S'. 
However, the argument used in section 2.2 for /9g 2 (i) is also valid for R(t',t), and implies 
that we can obtain R(t', t) by using the unsubtracted heavy-quark propagators. In order 
to compute R(t', t) we have used single and double cubic smeared interpolating operators 
J = hjo^q, with smearing size L s = 7, by using the heavy and light quark propagators 
rotated into the Coulomb gauge. L s = 7 was found to be the optimal value of L s for 



isolating the lightest meson state at (3 = 6.0 25, 26]. 

The procedure to extract operator matrix elements is standard. It is the same as the 
second method that we used in the previous subsection to determine C2- At fixed t, we 
study the behaviour of the ratio R(t',t) as a function of t' , searching for a plateau in t'. 
^bare ^ g dgf^gd D y weighted average of the data points in the central plateau region, 
if this exists. We will take as our best determination of Aj" 116 , the value evaluated in 
a time interval where the ratio R(t',t) appears to be independent of both t and t' . In 
addition, we have to require that the lightest state has been isolated. With the smeared 
sources used in the present case, we know that this happens at a time distance (t — t')/a 
and t'/a > 4 — 5 from the source. This implies that the total time distance t/a for 
R(t',t) has to be at least 8-10. Moreover, using (t — t')/a and t'/a > 4 — 5, we eliminate 



Figure 6: The ratio R(t',t) at t/a = 8 as a function of the time t' . We show the value of 
^bare _ j^f / a =3 ,5 R{t' , t) (full orizontal line) and the relative band of error (dashed orizontal 
lines). 

the contact terms, which on the basis of the discussion in the previous subsection, cf. 
fig. H are expected to be present up to distances of order 2-3. As an example of our 
results, we show in fig. ^ the ratio R(t',t), at t/a = 8, as a function of t', the time 
at which the kinetic operator is inserted. With the present statistical errors, it is not 
easy to identify the plateau region^. If we assume that we can use the central points 
(t' = 3, 4, 5) to extract the value of the matrix element, we obtain for the unrenormalised 
value a 2 X^ ie = -0.72 ± 0.14. This implies that there is a large numerical cancellation in 
the subtracted kinetic energy, a 2 X\ = a 2 Aj )are — C2 = 0.1 ± 0.14, cf. eq. (p4j). Due to the 
large statistical and systematic errors and to the difficulty in the clear identification of 
the plateau, it is not possible to obtain a value for the renormalised kinetic energy from 

At larger time distances, t/a > 10 the errors are even larger. 



this simulation. We can only impose the loose upper bound Ai < 1 GeV 2 . Notice that 
in order to reduce the statistical error to 0.1 GeV 2 , we need a sample about 50-100 times 
larger than our current one, corresponding to 1500-3000 gluon configurations. Moreover, 
we would eventually also like to be able to extrapolate the results to the chiral limit. For 
these reason, we are implementing the method described in this paper on the 24 Gigaflops 
APE100 computer. 

One could argue that the subtraction is not really necessary, since the effective theory 
on the lattice does not have renormalons. Even though this is indeed true, the difficulty 
in the determination of corrections of order I/tuq related to the kinetic energy operator 
would remain the same. The argument goes as follows. The bare kinetic energy operator 
has a very large matrix element a~ 2 x (a 2 A] )are ) ~ 2 2 x (—0.72) GeV 2 = —2.88 GeV 2 , while 
one expects a correction due to the kinetic energy of the heavy quark of the order of the 
squared Fermi momentum p 2 F ~ Aq CD ~ 0.1 GeV 2 . Thus the huge contribution of the 
matrix element of the bare operator has to be compensated by the corresponding term 
in the coefficient function of hh. This require an extreme accuracy in the perturbative 
calculation of the coefficient function. This remains true in the subtracted as well as in 
the unsubtr acted case. 



4 Conclusions 

In this study we have shown that the method for the non-perturbative renormalisation 
of the lattice operators h(x) D± h(x) and h(x) D 2 h(x), proposed in ref. [18|, is feasible in 
current computer simulations. We have been able to obtain the subtraction constants of 
the operators h(x) D 2 h(x) and h(x) h(x) with a small statistical error (in the former 
case, particularly for the constant C2 which is needed for many physical applications). 
The binding energy of the B-meson, A, has been also calculated with an error of about 
15% and was found to be significantly smaller than other estimates, based on different 
definitions pl| , p2[ . We have also computed the kinetic energy of the heavy quark in the 
B-meson. With our current statistical sample we can only impose the bound on Ai < 1 
GeV 2 , on the matrix element of the kinetic energy operator. This is due to the large 
numerical cancellation when the counter-term is subtracted. We are planning to improve 



the precision of our results by using a much larger sample of gluon configurations, and 
hopefully to obtain a significant result for Ax- 

Our preferred determination of A was based on the behaviour of the heavy quark 
propagator at large times. It is important to verify the validity of the condition (||) by 
extending the calculation of 6m(t) to larger values of t. This requires a high-statistics 
simulation on a large lattice, and we are currently undertaking such a study. The results 
will be reported elsewhere. 

The present study concerned some important matrix elements which appear in the 
HQET. We were able to determine A, defined in different prescriptions, with good preci- 
sion. This encourages us to extend the calculation to other matrix elements which appear 
at 0(1/toq), and beyond, in the HQET. The main limitation to the matching to the full 
theory is due to the fact that the relevant Wilson coefficients have only been computed 
at first order in a s . We are planning to extend these calculations to higher orders. One 
may also extend the present approach to matrix elements which appear at higher orders 
in other important operator expansions, such as the non-leading twist operators in deep 
inelastic scattering or higher dimensional condensates used in QCD sum-rules. 
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